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The calculation of finite temperature ESR spectra for concrete specified crystal configurations is a very im¬ 
portant issue in the study of quantum spin systems. Although direct evaluation of the Kubo formula by means 
of numerical diagonalization yields exact results, memory and CPU-time restrictions limit the applicability of 
this approach to small system sizes. Methods based on the time evolution of a single pure quantum state can 
be used to study larger systems. One such method exploits the property that the expectation value of the auto¬ 
correlation function obtained for a few samples of so-called thermal typical states yields a good estimate of the 
thermal equilibrium value. In this paper, we propose a new method based on a Wiener-Khinchin-like theorem 
for quantum system. By comparison with exact diagonalization results, it is shown that both methods yield 
correct results. As the Wiener-Khinchin-based method involves sampling over thermal typical states, we study 
the statistical properties of the sampling distribution. Effects due to finite observation time are investigated and 
found to be different for the two methods but it is also found that for both methods, the effects vanish as the 
system size increases. We present ESR spectra of the one-dimensional XXZ Heisenberg chain of up to 26 spins 
show that double peak structure due to the anisotropy is a robust feature of these spectra. 
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I. INTRODUCTION 


Quantum spin systems have attracted interests for decades because they exhibit various nontrivial behavior due to quantum 
mechanical effects. In particular, in low dimensions, quantum fluctuations due to non-commutativity of the spin operators and/or 
competition among the interaction (frustration) play an important role and various novel concepts, such as the valence-bond solid, 
resonating-valence bonds, and magnon Bose-Einstein condensation, etc. have been developed. One important topic is the effect 
of nonmagnetic defects. In a spin 5=1/2 antiferromagnetic Heisenberg chain, quantum fluctuations prevent the spins from 
being ordered, even at T = OK, and the ground state is nonmagnetic. However, nonmagnetic defects break the translational 
symmetry and polarize the surrounding spins cni . Then the one-dimensional system is described by an open-ended spin chain. 
An example of such a system is the Pd doped chain Sr2CuO3|0-01- 

ESR is one of the major tools to study the effects of defects in spin systems. Modeling the ESR spectra of intrinsic defects in 
spin chains is an important problem for which data for finite but rather long chains are necessary. In particular, the parameter 
dependence of concrete ESR spectrum for a specified system is of great interest and the temperature dependence of the ESR 
spectrum provides a lot of information about the spin ordering. To study these aspects theoretically, the explicit form of inter¬ 
actions and spatial configuration of magnetic ions in the lattice play an important role and it is necessary to study microscopic 
models, that is we should calculate the ESR spectrum for specific quantum spin Hamiltonian. The most direct way is to calculate 
the Kubo formula HU] by making use of the eigenvalues and their eigenvectors obtained by diagonalization of the Hamiltonian. 
The first attempt has been made to study the Nagata-Tazuke phenomena in a one-dimensional Heisenberg chain of 8 spins 
with dipole-dipole interactions d. In this work, the dependence of the spectrum on the angle between the static field and lattice 
direction was reproduced. Moreover, antiferromagnetic resonance d and also the appearance of a resonance forbidden by the 
Dzyaloshinskii-Moriya interacti onlfOl were studied. The structure of the ESR spectrum of a spin ring at high temperatures 
was investigated up to 16 spins lU^ . Complementary, Oshikawa and Affleck have developed a held theoretical approach by 
making use of the exact autocorrelation of Heisenberg chain, and have successfully analyzed low temperature properties for long 
(infinite) chains S. 

Obviously, the application of this method of exact diagonalization (ED) is limited to small systems for which we can obtain 
all the eigenvalues and their eigenvectors. Eor a system of A 5 = 1 /2 spins, we need the memory of order 2^^. By making use 
of the symmetries of the system we may reduce the dimensions of the block diagonalized Hamiltonian. Even with these efforts, 
in practice, the system size is still limited to about 20 spins. 

The restriction imposed by the exact diagonalization approach can be alleviated by computing the autocorrelation function 
(AC) from the time evolution of a pure state .The ESR spectrum is obtained by Eourier transformation of the autocorre¬ 

lation function of the transverse magnetization. In the AC method, a pure state evolves in time according to the time-evolution 
operator the action of which can be computed efficiently by means the Chebyshev polynomial lllSl - l^ or the Suzuki- 

Trotter-product-formula method ifTsl |2^ . As the memory required to store one pure state or all eigenstates is of the order 2^ and 
2^^, respectively, the time evolution method allows us to study systems that are twice the size of those that can be studied by 
exact diagonalization. Eor compute thermal averages we need, in principle, to take the average over the initial states. However, 
it is known that the expectation value of a quantity A for a single random state |<I>) yields an estimate of the trace of A, that 
is (<I>|A|<I>) ~ Tr A, which becomes more accurate as the size of Hilbert space increases O^ . Because of this fact, there is no 
need to compute traces of matrices to obtain the thermal equilibrium average dS I 22 I 424 I] . This approach has been used to 
study the temperature dependence of the total amplitude of ESR spectrum for the single molecular magnet V 15, consisting of 15 
5=1/2 spins m. Although this method to compute the ESR spectrum has considerable potential, properties of its applications 
to large systems have not yet been scrutinized in sufficient detail. An important issue, which we address in this paper, is that 
because the time-evolution method necessarily yields data for a finite time interval only, it is important to study how the ESR 
spectrum, that is the Eourier transform of this data, depends on the size of the time interval. 

In the AC method, the spectrum is obtained from the autocorrelation function ))eq. In ESR experiments, one measures 

the time evolution of the magnetization (t). The relation between the AC and the spectral density (f) is given by the Wiener- 

Khinchin (WK) theorem. This theorem relates the spectral density of the dynamics of a quantity and the Eourier transform of the 
autocorrelation function of the same quantity. In the present paper, we propose a method to directly compute the ESR spectrum 
from the time evolution of M^{t) by exploiting the idea of WK theorem. In quantum systems, however, the definition of the 
magnetization dynamics is somewhat tricky and we therefore develop a quantum version of Wiener-Khinchin relation, i.e., an 
explicit relation between Eourier transform of the autocorrelation function and the spectrum density in the quantum case. We 
call the approach based of this idea the Wiener-Khinchin (WK) method. Because the thermal average of the magnetization in 
the transverse direction (M^(f ))eq is zero, and the ESR signal is proportional to the average of square of the Eourier transform of 
over many realizations of the random state |<I>) ,the reasoning about the convergence of this average as a function 
of the dimension of the Hilbert space does not hold. Therefore, we present a detailed analysis of the distribution of the sampled 
data, and show that the distribution converges with a finite variance, independent of the size of Hilbert space. Because of the 
latter property, also for a large system, it is necessary to perform ensemble averaging of the data which renders the computational 
efficiency of the WK method less than the one of the AC method. However, the WK method gives additional information about 
effects of flniteness of time domain. In particular, the Gibbs oscillations in the WK method are positive only whereas they can 
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be negative in the AC method. We find that as the size of system increases, for both methods, the effects of finite time interval 
disappear and the ESR spectra are very similar. For large systems, the AC method is the most efficient one because a single 
sample of a random state yields accurate result. We use both methods for applications to large systems which have not yet been 
studied in previous works. In particular, we compute spectra for one-dimensional XXZ model up to 26 spins and show that 
the double peak which was found in earlier work lo is almost independent of the system size, indicating that the double-peak 
structure will be present in the thermodynamic limit. 

The outline of this paper is as follows. A brief overview on the methods previously used is given in Sec [III In Sec. [Till we 
introduce the new method motivated by the Wiener-Khinchin theorem, and study statistical properties of the method in detail 
in Sec, mid In Sec. lIVl we show the distribution of sampled data of the spectrum obtained by the AC and WK methods. The 
total amplitude of the spectrum, which is fundamental information in ESR studies, is discussed in Sec.[V] Section [VTl presents 
applications of the methods to large systems. Summary and discussion of related problems are given in Sec. lVIIl Appendix A 
discusses the effect of the finite time interval on the spectrum for the both methods and Appendix B gives a detailed analysis of 
the statistical properties of the thermal typical state for the both methods. 


II. RECAPITULATION OF EXISTING METHODS 
A. Kubo formula and ESR spectrum 

The ESR spectrum is given by the Kubo formula|0[Ml- The imaginary part of the dynamical susceptibility z"(®) read^ 

Z"(®) = ^(1 1“ (M"(0)M"(0)eqe-“'df, (1) 

and the ESR absorption spectrum is given by 

I^ico) = -^x'\oy), ( 2 ) 

where Aq is the amplitude of the external field, and we adopt 

N 

M^t) = ^ S-, (3) 

1=1 

(•)eq = Tr[ • (4) 

For the numerical calculation the ESR absorption spectrum, several methods have been developed. There are essentially two 
types of methods: (1) Exact diagonalization method iflolfT^ and (2) Time-evolution of the autocorrelation function method ifl^ . 


B. Exact diagonalization (ED) method 

The most direct calculation of the Kubo formula uses the set of eigenvalues and eigenvectors {£„, |n)}^=i obtained by solving 
the eigenvalue problem 


J^f\n)=E„\n) (5) 

for the hamiltonian D denoting the dimension of the Hilbert space. In practice, we solve Eq. @ by numerical diagonaliza¬ 
tion. The autocorrelation function (M^(0)M^(f))eq is expressed as 

(M^(0)M^(f))eq = ^ /Z, (6) 


where Z is the partition function 



( 7 ) 


* Throughout the paper we denote the temperature by j3 ’in order to avoid confusion with the end of the time interval T. 
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The Fourier transform of Eq. (|6l) reads 

r (M-’^(0)M"(f))eqe-‘“'dr = £ | {m\M-'‘\n)\h-^^'-2n5 (co - {E^ - E„)) jZ, (8) 

Jm.n 

where we used the definition 

5{(0-{E^-En)) = ^j_ (9) 

The imaginary part of the dynamical susceptibility reads 

= (to(10 ) 

m,n 

where 

= (0,n,n=E„,-E„. (11) 

It is sufficient to consider (Om,n > 0 only since we are interested in the absorption, not the emission. Note that x"{^) > 0 for 
0 ) > 0 . 

In the exact diagonalization approach, the spectrum consists of a finite sum of delta functions. Therefore, to draw the spectrum, 
we have to represent each delta function by a regular function of certain width. For example, we may replace the delta function by 
a Gaussian or simply use a histogram representation. The results of the exact diagonalization are exact, usually close to machine 
precision, but we need to know all the eigenvalues and corresponding eigenstates. Therefore, we need memory of the order of 
for a matrix of the size D = 2^ for systems of N S = 1 /2 spins. If the system has symmetry, we may reduce the size. For example 
if the system conserves M"’, then D is reduced to nQv/ 2 +m:- Moreover for ESR spectra, only the uniform mode is relevant and 
therefore only the fully symmetrized states are necessary, which also allows a reduction of D. However, the memory limitation 
prevents us from studying more that N = 20 S = I/2 spins. As an illustration, we use the exact diagonalization approach to study 
the one-dimensional spin-1/2 XXZ model in a static magnetic field H and compute the response of the magnetization along the 
X axis to the AC-field A(f). The hamiltonian of the system is given by 

N N N 

r=l !=1 !=1 

where we impose the periodic boundary condition a = x,y,z unless otherwise mentioned. Hereafter we adopt Kelvin 

as the unit of energy and we set gjig = 1. 

In Fig. [T] we show an example of spectrum obtained by exact diagonalization method for an antiferromagnetic Heisenberg 
chain with A = 16 with the static field El = 5K at the temperature = lOOK. The notation and the parameters of the model 
are set to be the same as those in the model studied in the previous studv lfT^ {Jx = Jy= IKjT^ = 0.92K). 

The spectrum consisting of an ensemble of the delta peaks (Eq. ( fTOl i) is plotted as a histogram with small bins Aw. In the 
exact diagonalization (ED) method, we are free to choose any value of Aw. Anticipating for the comparison with the two other 
methods, we take Aw = 2nlT where T = df x df = 0.5 and = 16384 is the number of data points. 



FIG. 1. Spectrum for a ring of A = 16 spins. The histogram is obtained hy the exact diagonalization method. 
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C. Autocorrelation function method 

According to the Kubo formula (Eq.([Tli), the spectrum is given by Fourier transform of the autocorrelation function. Because 
of the symmetry 




(13) 


where * denotes complex conjugate, we have 




e-‘“'df = 2Re 




e-'“df 


(14) 


where Re[-] denotes the real part. Therefore, to compute Eq. (fT4l) . it suffices to have AC data for the interval [0, T]. 
We obtain the spectrum by discrete Fourier transform (DFT) of the AC data 

jfiT 

fi^m) = ^ (^m))eq; Un — j HI = 0, 1,2, * • • ^flt 1, 

nt 

where «r is the number of data items. The DFT of Eq. (fl4l i is given by 

j- «f— 1 


(15) 


n. 


nt—\ / mT\ irlr 

E / — , ^=-„„-„, + l,...,0,l,2,...,u,-l. 

n \nt J T 


(16) 


Note that the absorption and emission spectrum correspond to (Ok > 0 and cOk < 0, respectively. The maximum angular frequency 
{2n times Nyquist frequency) that can be represented by the DFT Eq. (fThl l is given by niit/T. For a magnetic field H ~ 5K, the 
main contribution to the absorption spectrum peak comes from a narrow peak around co ~ 5K. In our numerical work we take 
df = r jut = 0.5. Therefore the largest value of omega {2k times the Nyquist frequency) is {2k/T) x {ntl2) = n/dt = 2 k m 
6.28 > 5, which is large enough to cover the full spectral range. 

In Fig. |3left), we show the spectrum for a system of A = 6 spins at the temperature = lOOK, obtained by using a 
time series of iif = 16384 items. Clearly, the spectrum is suffering from fine oscillations with negative values values, which 
is called the Gibbs oscillation ll^ . This oscillation is due to the fact that the range of the time integral is finite and that there 
are eigenvalues that do not not exactly match one of the cOk- This artifact can be suppressed by employing a suitable window 
function lE^ . for instance a Gaussian window (see Appendix A). As shown in Fig.|2](Right), the Gibbs oscillation is reduced and 
the result is consistent with the one derived with the exact diagonalization method with the resolution of 0{a) {a is an artificial 
parameter which determines the resolution of the spectrum. See Appendix A for the precise definition). 




FIG. 2. The spectrum obtained by exact diagonalization method and by the time-evolution method: N = 6, temperature = 100K,Ff = 5K, 
sampling time dt = 0.5, the number of data items n, = 16384, and the to-mesh njT = K/{nt x df) (blue line). The exact diagonalization result 
(green line), uses an co-mesh Aco = 2;r/16384, which corresponds to the frequency n/T used in time-evolution method. Left: the spectrum 
obtained by exact diagonalization (green line) and by time-evolution method without window functions (blue line). Strong Gibbs oscillation 
is found. Right: Same as left except that a Gaussian window function with standard deviation a/T = 1 /9>\92 = 0.00085 (a = 7) is used. 
Clearly, the Gibbs oscillation is suppressed. 
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1. Thermal typical state method 

In theory, the calculation of the thermal equilibrium expectation value of involves a trace operation. This requires 

of the order = 2^^ operations, which becomes prohibitively large if N is of the order of 20 or larger. Fortunately, for larger 
systems, we can obtain accurate estimates of thermal equilibrium averages by making use of the so-called thermal typical state 
ifl^fl^l^l^ . also called “a Boltzmann-weighted random vector” |[l^[l3] or “a canonical thermal pure 
quantum state” ll^ , where |<I>) denotes a random state on the D-dimensional hypershere. In essence, we have 

TrXe^l^-^ 

= = (17) 

We briefly explain the ideas behind this method. Let be an arbitrary set of complete orthonormal states of the Hilbert 

space of the system. Using the complex-valued random variables {^n}n=v introduce a random vector: 

l‘J>)= (18) 

n=\ 

Note that this construction is independent of the choice of the basis set {|n)}^=i. 

Unlike in [A. Hams and H. De Raedt, PRE (2000)] , to simplify the mathematics, we will work with independent complex¬ 
valued Gaussian random variable for each index n. They have identical Gaussian distributions with mean zero and variance 2(7^ 
for all real and imaginary parts of the variables. By E[-], we denote the expectation with respect to the multivariate Gaussian 
probability distribution: 


D 


D 


P{^\ , • • • , n ^a) = Y[ 

a=\ a=\ 


2na^ 


d(Re^„)d(Im ^„). 


(19) 


The following relations hold: 


E[^„]=E[^J=E[^„4]=0, 

E[^:^=2a^8a,p, 

+E[4:^JE[4*^,] = 4a" , 

which are used in the following calculations. 

For any matrix X we have 

E[(<I>|A|4>)] = f E[^:^,,]{a\X\p)=f^E[^:^^,]{a\X\a)=2a^TrX, 

a,p=l a 

and because ('1>|A"|<I>) = (OlAlO)*, the corresponding variance is given by 
Var((<I>|A|<I>)) = E[|(4>|A|4>)p] - |E[(4>|A|<I>)] ^ 

= f E[^;^^4^;](a|A|p)(h|A|^)*-4a"|TrAp 

a,p,b.q=l 

= 4a" Y. ({a\X\a){b\X\b)* + {a\X\b){a\X\b)*^-4a^\TrX\^ 

a.b=l V / 

= 4a"TrAAl. 


The relative standard deviation is dehned by 


RSD2((<I>|A|<I>)) 


Var(((I>|A|(I>)) 

MwW 


TrAA" 

\TrXf 


TrAA" 

(TrA)(TrAt)' 


Next, we introduce the thermal typical state |<I>| 3 ) for the system at an inverse temperature j3: 


|4>^)=e 


( 20 ) 


( 21 ) 


( 22 ) 


(23) 


(24) 
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where (<I>| 3 |<I>j 3 ) = (<t>|e is an approximation to the partition function Z^. As shown in Appendix B, the RSD of 

(<I>|e^^‘^|4>) vanishes with D as 

RSD ((Op I Op)) = RSD((0|e-^‘^|0)) < (25) 

Thus, hereafter we only estimate the RSD of the numerator for the thermal average. In our calculations, we take averages of the 
dominatorand the numerator with respect to the same samples {|Op)}: 




(26) 


where N is the number of samples. But if we take the average of the ratios of samples, we obtain the same result in the limit 
D (see Appendix B): 


/Y\ ~ 1 V 


(27) 


The autocorrelation function (MW*(f))eq is given by 


(M"M"(f))eq = 


E[(Op |Op)] 

E[(<J>p|Op)] 


(28) 


The RSD of the autocorrelation function also vanishes with increasing D (see Appendix B). Thus, for a sufficiently large system, 
a few samples suffice to estimate the autocorrelation function. Here it should be noted that strictly speaking Eq. (fOl) does not 
hold for the quantity (<t>p |M^e''^^M''e^'‘^'|<I>p), but for sufficiently large D, Eq. (|28]) justifies using Eq. (fOl) . 

Eor the autocorrelation function {M^M^{t))eq, we need to construct the states 


|A) = e-'-^'|4>p), and |B) = e-‘-^'M^|<I>p), 


(29) 


and then obtain the expectation value of the autocorrelation as 

(Op iMV-^^M^e-'-^'lOp) = (B|M^|A). 


(30) 


To perform the operations e P-^/^ and e on the states, we may use the Chebyshev method or the Suzuki-Trotter decom¬ 
position method. In the present study we used the Chebyshev method for e^^'^/^jO) and mostly used the product method for 
e-i^'lOp). 


III. WIENER-KHINCHIN METHOD 


In this section, we propose a new method that makes use of the Wiener-Khinchin theorem. In the autocorrelation method 
(see previous section) the spectrum is estimated from the autocorrelation function {M^M^(t))eq- In experiments, the spectrum 
is obtained from the record of time evolution of magnetization The relation between the autocorrelation and the Eourier 

transform of time evolution of a quantity 2f (f) is given by the Wiener-Khinchin relation, that is, the relation between the spectrum 
of the fluctuation in time of a quantity X (f) and the Eourier transform of the autocorrelation function {X (O)^(f)). 

We explore the possibility to obtain the spectrum from the dynamics of the quantity itself, i.e. (Op \X (f) |Op) from an initial 
state |Op). In this method, we use the time-evolution of the state as in the previous subsection and therefore we can study large 
systems as well. 


A. Wiener-Khinchin theorem 

Eirst, we briefly review the Wiener-Khinchin theorem. Eor a time-dependent quantity X{t), we define the autocorrelation 
function R{t) and the spectral density S{co) as 

R{t) = {X{0)X{t))= lim i rx{z)X{t + z)dz, 

1 Jo 


(31) 
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and 


5(®)=limE^, X^(®)= rx(r)e-™df, (32) 

T^oo I Jq 

respectively. The Wiener-Khinchin theorem tells us that the Fourier transform of the autocorrelation function equals to spectral 
density: 

G(®)= [ R{t)e-^“‘dt = S{co). (33) 

Note that we assume that the process ^(f) is stationarv i^ . 


B. Dynamics of the magnetization 


Next, we apply the idea of Wiener-Khinchin relation to the Kubo formula. By making use of the relation (l3?t we obtain 
G{co) from the spectral density S{co). Here it should be noted that M^{t) is an operator and the definition of {M^{t)) is tricky. 
Definitely, (M^(f))eq is time-independent, in fact zero in the present case, and therefore we cannot extract S{(£)). 

In the following, we propose a method to obtain G{(£)) from {M^{t)) in a quantum mechanical system. This approach is 
motivated by the Wiener-Khinchin theorem but we do not use the relation (l33l l directly. In time-domain methods, we use the 
notion of the thermal typical state (see above) and we use this notion here once more to derive a formula to obtain S{co) from 
time evolution of {M^{t)) ll2^ . 

First we prepare a thermal typical state {‘Pp) as an initial state. The expectation value with this state gives the thermal average 
at the inverse temperature j3. For (M^(f))eq = 0 and the expectation value E[(<I>j 3 |M'’^(f)|<I>| 3 )] is zero, but, in general, 

(<I>| 3 |M^(f)|<I>j 3 ) is not zero. When we calculate the time evolution of the sampled state we can extract information for the 
spectral density 

Mj{(0) = £ {^p\M'^{t)\<i>p)c-''^^dt. (34) 

Although E[{<i>p\M^{t)\<Pp)] is zero, the expectation value of the spectral density E[|M^(a))p] is not zero and, as we show 
below, we can obtain the ESR spectrum from the latter quantity. 

Taking the basis to be the eigenstates of the system, {<Pp \M^{t)\<i>p) is expressed in terns of the random numbers as 

(4>^|M"(f)|<I>^) = (35) 

m,n 

where \n) and £„ are the eigenvector and its eigenenergy of the system Hamiltonian (Eq. Q). The Eourier transform Mp {(o) 
(Eq.(l34]i) is expressed as 

Mj{(0) = {(0- {E,n-En)) {m\M^\n), (36) 

m,n 

where we introduced the function 

sin — 

(37) 

71(0 

Note that we cannot extend the range of integration as in the AC method (Eq. (O) because for a given \‘i>p), {‘i>p |M^(f)|<I>j 3 )* = 
p\M'^{—t)\<i>p) does not hold in general. Eor the same number of time steps, this implies that the frequency resolution 
Aco = 2%lT is twice as larger as for the AC method. 

As a next step, we calculate the average over the random initial states of the quantity 






mMm'.n' 


X 4;r^5^ {(O - {Em - En)) 5^ {co - (£„/-£„/)) {m\M^\n){n\M-'‘\m). 
Using Eq. ( l20b and the fact that («|M^|n)=0, we have 

E[\Mj (a))|2] = 2 a^Te-P‘^'£e-^P^'’ 27 t 5 ^ {co - (£„-£„)) |(m|M"|n)|2, 

m,n 


(38) 


(39) 
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where we used the relation 


{8'^{(0-{E,n-En))f-^5'^ {(0-{E„-En)). 


The spectral density reads 




E[|M^(®)|2]/z2 

lim- - - - 

T —>oo T 

^ c-^P^''27z5 (o) - (£„, - £„)) | (m|M"|n) | Vz^, 
m,n 


where 


(40) 


(41) 

(42) 


Z^ =E[(<I>;3|‘I>/3)]. 

The Eourier transform of the autocorrelation function is given by 

Gp{(0) = '£e-P'^"2n8 {(0 - (£„ -£„)) | {m\M^\n)\yZp. 

m.n 

By comparing (1421) and (l44l) . we obtain a Wiener-Khinchin-like relation for the transverse magnetization: 

’jl 

Gp{(0) = 2a^-^c^i:p^2{(0). 

The imaginary part of the dynamical susceptibility of interest is given by 

X"{0)) = (T^^sinh ^/2i(o}- 

Here it should be noted that we need to calculate the quantities of j3/2 (not j3) to obtain the ESR spectrum of j3. 

In Eig.[2 we show a comparison of the spectrum obtained by the exact diagonalization method and by the WK method. 




(43) 


(44) 


(45) 


(46) 


FIG. 3. Absorption spectrum obtained by the ED method (green line) and by the WK method (red line): N = 6, = lOOK, H = 5K, and 

the mesh of the frequency A® = 2;r/8192. In the ED method (green line), we take the mesh of the frequency A® = 2;r/8192, which is the 
interval of the frequency l^jT in the WK method. (Left) The Wiener-Khinchin spectrum is obtained by calculating Eq. l l39b exactly. (Right) 
The WK spectrum is obtained by averaging 10000 samples of thermal typical initial states. 


Erom Eig.[3it follows that the spectrum is well reproduced without any window procedure, in contrast to Eig.|2](Left) obtained 
by the AC method. This fact is attributed to that the finite time effect (Gibbs oscillation) in the WK method gives positive and 
smaller artificial peaks in the spectrum (see the Appendix A). 
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C. Sample distribution of data 

In principle, methods that are based on the thermal typical state require sampling over the random states to obtain data for 
equilibrium. However, as mentioned above, the data obtained by the thermal typical state converge to the thermal equilibrium 
value as the size of system increases exponentially. That is, the deviation decreases proportionally to where D is the 

dimension of the Hilbert space. In AC method, the data are obtained as expectation values of and for large system the 

convergence to the equilibrium value is so fast that one initial state suffices to obtain accurate results (see Appendix B for the 
details). On the other hand, the WK method extracts information from the quantity Mp (cu) whose expectation value is zero, and 
therefore we have to investigate the fluctuations of the squared variable \Mp (w) p because it is not obvious that these fluctuations 
will decrease rapidly as D increases. 

In this section, we study properties of sample distribution of the data. To characterize the distribution, we study the variance 
of \Mj{co)\^: 


Mp{(o) = {<P\A{(o)\<i>), 

1 


To this end, we introduce 

where 

and 


n 2 




e-'“'ei-^'M'’^e-i-^'df. 


r /2 


(47) 

(48) 

(49) 

(50) 


which is proportional to Eq. d^ . i.e., Mp{co) = Mp{(o) x The reason for shifting the interval from [0,T] to 

[—T /2, T /2] is to assure that Y (o) = Y^(—co) and 

A,^^A{co)ap = {A{-co),,r = {AHco)pay = ( 51 ) 

where we choose the states \a) to be the eigenstates of M’. With this choice, {a\M^\a) = 0 and Aaa = 0 which simplifies the 
calculations significantly. Using these quantities. 


a,b,p,q 

and, using the identities Eq. ( |20] |. the expectation of Eq. ( |52] | is expressed as 

E[|M^(a))p] = = Td'^TrAAt. 

a.b 


(52) 


(53) 


Hereafter, for simplicity of notation but without loss of generality, we choose 2a^ = 1. Note that in Eq. (l5^ . we did not yet 
include other terms in ( |46] |, i.e., Zp. But, as we mentioned above, the variances of them are small and thus we only concern the 
variance gTli. 

We calculate the first term in Eq. (l47l i by making extensive use of the properties of Gaussian random variables. We have 

a.b.c.d p,q,r,s 


+E[^;^,]E] +e[4:4]e(54) 


a,b,c,d p,q,r,^ 


The first term in the curly brackets vanishes because E[^*^p] = and Aaa = 0. Interchanging the summation indices {r,q) 
and {s,q) in the third and fourth term, respectively we obtain 


= TTrAAAU'l' + jTrA^j^ + 2 (TrAA^)'" +2TrAAUA^ 


(55) 
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and using Eq. (l53T l we finally obtain for the variances 


Vai(\Mj{a))\^^ =4TrAAAU'*' + \TrA^\^ + {TrAA'^f + 2TrAA'^AA'^. 


(56) 


We want to find bounds to the estimate of the variance Eq. ( l56l l. Let us denote the (non-negative) eigenvalues of AA^ (A^A) 
by x| (yl = x\), that is AA^|x^) = x||x^) (A^Ajjit) = yl\yk))- We assume that the eigenvectors of AA^ and A^A are normalized. 
Then we have 


>^x^ = TrAAUAt, (57) 

TrAAUU = X^xfyf |(x,|y,)|2 < ^xfy^ = (XrAA^) (XrAU) = (TrAA^)^ (58) 

k,l k.l 


As {X,Y) = TrX^y defines a scalar product, by the Schwarz inequality |(A’*',A)p = |TrA^|^ < |TrAA’^|^. Noting that the last 
term in Eq. (|56] | cannot be negative, putting all this together we find 



1 < 1 + 


TTrAAA'^A^' + jXrA^ |^ + 2XrAAUAl' 
(XrAAt)^ 


< 8 . 


(59) 


Thus, we find the variance is bounded as 


Var 

1 / 


E 


2 


(60) 


indicating that in the WK method, it is necessary to average over different initial states, also for large systems. This is in sharp 
contrast to the case of the AC method in which the variance vanishes exponentially with the system size. But Eq. (l60l) also 
indicates that the variance is bounded, independent of the system size, and therefore we can obtain good estimates for the mean 
from finite samples. 

We emphasize that the variance of the WK method does not depend on D, that is, the distribution shows a kind of typicality. 
Thus we can obtain the correct mean with aimed precision by sampling, regardless of the system size. 


IV. DISTRIBUTION OF AVERAGE VALUE AND CONVERGENCE IN DISTRIBUTION 


In this section we investigate the sample distribution of x(ti)) data. Eor concreteness, we choose the value of co to be 4.9916K 
at which the spectrum has a maximum (peak). 

Eirst we show the distribution of data in the AC method for the system A = 10 where D = 1024 and D ^ Ri 0.001 is rather 
small. The histogram of this data is given in Eig. |4|Left). The exact value is denoted by the black line. The mean of the 
distribution is shown by the bold blue line. The standard deviation (square-root of the variance) is given by the arrow. Clearly 
Eig.llfLeft) shows that the data is distributed around the correct value with RSD=0.31025K. As we already menioned several 
times, as the size of the system increases the variance will vanish as D *. 
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FIG. 4. (Left) The distribution of values over 1000 samples obtained by the AC method at to = 4.9916K for N = 10. The mean (=4.8508K) is 
given by the blue line, and the standard deviation (=1.5049K) is shown by the arrow. The exact value is denoted by the black solid line whose 
CO = 4.4666K. (Right) The distribution of values over 1000 samples obtained by the WK method at ® = 4.9916K for N = 10. The mean 
(=4.01579K) is given by the blue line, and the standard deviation (=3.9749K) is shown by the arrow. The exact value is denoted by the black 
solid line whose co = 4.4666K. 

The histogram obtained by the WK method is depicted in Fig.|4jRight). In contrast to the AC method, the distribution is of 
the exponential type. From the figure, we find the variance is about the same as the average, i.e., RSD=0.98982K. This value 
corresponds to the lower bound of the estimation dbOl l. This fact will be discussed in more detail below. 

We present the spectra for N = 10 obtained by the both methods compared with the exact one obtained by the ED method. 
Because we have the exact eigenstates and eigenvalues, we computed the results of AC and WK methods by employing Eqs. (l36l l 
and (l44li using the exact eigenstates and eigenvalues (not by sampling). 




FIG. 5. (Left) Spectra obtained by the AC method (blue line) by the WK method (red line) with nt = 4096. The column histogram is obtained 
by the ED method. The mesh of the frequency 2;r/4096 for the WK method and ;r/4096 for the AC method and the ED method. The line 
CO = 4.9916K is denoted by the black line. (Right) n, = 16384, and the mesh of the frequency 2;r/16384 for the WK method and ;t/16384 for 
the AC method and the ED method. 

Although both methods approximately reproduce the exact results, the AC method shows large fluctuations for both tit large 
and small. On the other hand, the WK method gives rather good agreement. 

In the case of the AC method, the obtained values have a narrow distribution around the mean, and it is known that the variance 
decreases as D In this sense, the AC method has an advantage. But the WK method also gives correct estimates by sampling 
and the obtained shape of the spectrum is close to the exact values. 


A. Special case 1: Energy gap selection 

In the previous section, we found that the relative variance is close to one, the minimum value of the bound Eq. dbOl l. In this 
subsection, we consider the reason of this fact. To derive Eq. dbOl) . we considered the general case. But we can also use the fact 
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that Aap (Eq. (ISTT i) is almost zero except when (oc^ — Ep. If we assume that there is a unique set {n,l,m) which satisfies the 

conditions = E„, — Ei = (O and {m\M^\n){m\M^\l) ^ 0 then all terms are negligibly small except for the term (Tr 

Therefore 

^Var(|Mp^(at)P) ^ (xrAAt)" 

E[|M^(a))|2] “(TrAAtf”^’ 

which is indeed what we have observed in our actual calculations. 


B. Special case 2: Isolated delta function 


We give the variance of the WK method for the case where the spectrum consists of well-separated delta functions which are 
given only by a single resonant pair of states, that is the non-resonance condition that if Em — En = — E„i then (m = m' and 

n = «') or (m = n and rd = «') is satisfied. Then, we can use the fact that approximately (x) = 5(x) and we have 

E[\Ml{co)\^]=2a^TQ-^^ {co-{Em-En))\{m\M^\n)\^. (62) 

The variance is estimated by making use of the relation: 

E = E AAn? ^ 4 (63) 


as 


E 


|M^(cu)r-E[|M^(a))pf 




(64) 


Thus, the relative variance is 3, a result which does not depend on D. 

In the AC method, the variance of sampled amplitude of the isolated delta functions is the same as the expectation value for 
the present case and RSD > 1. However, the isolated delta function in the spectrum is only relevant in small systems. For large 
system with large D, the amplitude of the isolated delta function is small and we can ignore it. 


V. THE TOTAL AMPLITUDE OF THE SPECTRUM 


The total amplitude of the spectrum, or the intensity, i.e. the integral over the absorption spectrum (Eq. ©), is given by 


/^ = 



(65) 


It is one of the fundamental pieces of information obtained in ESR experiments. Its temperature dependence has been calculated 
for the single molecular magnet Vis It^ . In this subsection, we demonstrate that P obtained by the WK method correctly 
reproduces the exact ED results. Details of the analysis are given in Appendix A. In the WK method, a delta function is 
represented by the sinc-function 




1 / sinwr/2\ 


K(0 


■ 


Because of the relation 


1 

T 



( 66 ) 


(67) 


the intensity P is obtained correctly by analytical integration. However, when the DFT yields the spectrum at discrete points co = 
Ikn/T, k = 0,l,2,-- - ,«? — 1 only. Thus, the integral over the co is given by a discrete sum. In Appendix A, we demonstrate 
that the discrete sum agrees with the analytical result. We also show that for the AC method in which the discrete points are 
given by 0) = kn/T, k = —tit,--- , —1,0,1,2,--- ,n, — 1 we recover the analytical result of the intensity, in spite of Gibbs 
oscillations. We calculate the results of AC and WK methods by evaluating the formulae for the expectation value by using the 
exact eigenvalues and eigenstate, as we did for Fig.|2]and[3^Left), i.e., without sampling over thermal typical states. In Fig.|6] 
we show that both the AC and WK methods reproduce the results of the ED method. 
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FIG. 6. Comparison of the intensities of the absorption as a function of the temperature obtained by the WK method (red cross) and by the AC 
method (blue open circle). The green solid curve is the exact value; N = 6, H = 5K. 


VI. APPLICATIONS TO LARGE SYSTEMS 


As an example of application of the methods to larger systems, we compute the spectrum by the AC and WK method for 
N = 16 with lit =16384. For this size, we can still compute all eigenvalues and eigenstates, hence we know the positions of 
the eigenfrequencies and the corresponding amplitudes. The ensemble of delta functions is turned into a histogram by using a 
frequency mesh A (0 = In/T, (T = 16384 x df (= 0.5)). 

In Fig. EtLeft), we present a comparison between the spectra obtained by ED and AC methods. The data obtained by the 
AC method shows sharp features because we did not use a window in the DFT procedure. As the statistical fluctuations on the 
autocorrelation are small, 10 samples suffice to find almost perfect agreement with the ED result. Therefore, we may conclude 
that as the size of the system increases, the effects of finiteness of the time interval reflected in e.g. Gibbs oscillations seems to 
be small whereas we found serious effects in the case A = 6,10. 



4.9 4.92 4.94 4.96 4.98 5 5.02 5.04 5.06 5.08 5.1 

co/7_. 



4.9 4.92 4.94 4.96 4.98 5 5.02 5.04 5.06 5.08 5.1 

co/7_. 


FIG. 7. (Left) Spectrum for N = 16. The histogram is obtained by exact diagonalization, points with error bars are the data obtained by the 
AC method with 10 samples and the time interval [0, T = 16384 x df],d/ = 0.5. (Right) Results of the WK method with 52 samples over the 
sme time interval 

In Fig.|7jRight), we show the comparison between the spectra obtained by the ED and WK method. In the WK method, we 
took 52 samples, and we did not use a window function for the DFT procedure. Here we again And a good agreement but the 
variance of ensemble average is about the same as the mean as we found in the case A = 10. We conclude that for A =16, the 
results of both methods are in good agreement those of ED. 

With above observation, we also obtained spectra for A > 20 by the AC method. Up to now, we adopted the periodic 
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boundary conditions (PBC). For a change and also because open boundary condition (OBC) is of interest because they occur 
in actual materials, we show results for N = 20 obtained for both boundary conditions. In Fig.[8jLeft) , the red curve (circles) 
is the spectrum obtained with PBC and the blue curve with triangles is that of OBC. We find that the spectrum of OBC has a 
shaper double peaks in the center, but the global shape is similar in both cases. Finally, we present spectra for N = 24,26 for 
OBC in Fig.0(Center) and (Right), respectively. We used n, = 16384 with df = 0.5 for = 20, and n, = 8192 with df = 0.5 
for = 24,26 and did not use any windowing procedure. From the size dependence, we conclude that the double peak structure 
will survive in the large size limit. 





4.9 4.92 4.94 4.96 4.98 5 5.02 5.04 5.06 5.08 5.1 4.9 4.92 4.94 4.96 4.98 5 5.02 5.04 5.06 5.08 5.1 4.9 4.92 4.94 4.96 4.98 5 5.02 5.04 5.06 5.06 5.1 

ai/J^ aUj. (B/y,. 


FIG. 8. Spectrum obtained by the AC method with the open boundary condition (blue line). (Left) N = 20. The red line denotes the spectrum 
with the periodic boundary condition. (Center) N = 24. (Right) N = 26. 


VII. SUMMARY AND DISCUSSION 


We have studied time-domain methods to compute the ESR spectrum through sampling of thermal typical states. Although it 
has been known that the autocorrelation function can be obtained with a few samples thanks to the typicality property, practical 
aspects of these calculations have not yet been explored. 

We proposed a method to obtain the spectrum from the motion of by making use of the Wiener-Khinchin-like relation. 
The dynamics of magnetization from sampled initial states provides an estimate the spectrum density. But in the quantum 
equilibrium state the expectations value of M^{t) and the sampled data for the spectrum density Mp{co) are zero. Then, we take 

thermal typical states as initial states ,and obtain the average of squared spectrum density \Mp (o) p, i.e., the Fourier transform of 
the autocorrelation function. To this end, we proposed a quantum mechanical version of the Wiener-Khinchin relation. However, 
it should be mentioned that the WK method does not directly relate to the experimental situation. Because (M^ (®))eq = 0, the 

averaged squared spectrum density |M^(a))p does not converge to a mean value but has a finite distribution, also for large 
systems. Therefore the convergence property of the thermal typical state does not apply in this case. We studied statistical 
properties of the distribution and found that form of the distribution converges in the sense of typicality. We derived bounds on 
the relative variance of the distribution. 

In time domain methods the effect of finite observation time is an important issue. We studied this aspect in detail for both the 
AC and WK method. It was found that for small systems, the AC method is suffering from severe effect of finiteness of T while 
in WK method the effect is practically suppressed. The WK and AC methods give complementary information for the effects of 
finite T and a comparison between them is useful to confirm that the obtained result is not affected by the finiteness observation 
time. We also found that for both the AC and WK method, the effect of finiteness of T decreases with increasing system size. 
As the system size increases, the efficiency of the AC method increases too. 

We presented the spectrum for a one-dimensional XXZ chain up with to A = 26 spins, where the spectrum shows a double 
peak which will be sustained in the thermodynamic limit. In our numerical work, we focussed on chains with an even number 
of spins. The ESR spectra of chains with an odd number of spins show behavior that is qualitatively different from the one of 
chains with an even number of spins, meriting a study in it own right. As the focus of the present paper was on time-domain 
methods rather than on specific applications, we relegate this study to future research. 

In summary, in this paper, we proposed the new method to compute the ESR spectrum by making use of mathematical relations 
discussed in Sec |III] We have studied one particular time-domain method based on the Wiener-Khinchin theorem but there are 
many other ways to obtain the spectrum density by making use of relations similar to Wiener-Khinchin theorem. Eor example, 
if we calculate E[{<i>p\M^\<Pp){<Pp\M^{t)\<i>p)], we can obtain a similar expression. Studying the relation among them is an 
interesting problem for future research. 
























16 


ACKNOWLEDGEMENTS 


The present work was supported by Grants-in-Aid for Scientific Research C (25400391) from MEXT of Japan, and the 
Elements Strategy Initiative Center for Magnetic Materials under the outsourcing project of MEXT. The numerical calculations 
were supported by the supercomputer center of ISSP of Tokyo University. We also acknowledge the JSPS Core-to-Core Program: 
Non-equilibrium dynamics of soft matter and information. 

Appendix A: Effect of flniteness of the observation time domain and window function 
1. Effect of Gibbs oscillation 

Because the observation time is finite [—7’,?’] and we have a finite number of points in time, the discrete Eourier transform 
(DPT) gives amplitude at the discrete points 


(Ok = Ao) X k, Aco = —, k: integer. 


(Al) 


On the other hand, the resonance points are not necessarily at those points. If the spectrum has a delta-function at co, this is 
expressed by several point at the discrete points in the spectrum obtained by the DPT s. We study this phenomenon for the 
AC and WK methods. 

Por a time sequence of data {fk}'k=-oo as the autocorrelation function, the spectrum of the process is given by the DPT: 



(A2) 



(A3) 


where A is the sampling interval. Por a set with a finite number, say 2n, the transformation reads 


P(a)) = A'^ /,e-‘“^^ = A £ /,/r,e-““ 


(A4) 



where 




otherwise. 


(A5) 


Expression Eq. (IA4t can be written as a convolution: 



(A6) 


where H{(o) is the DPT of {hk}k 



(A7) 


indicating that the spectrum we want is deformed by 



(A8) 


In Pig. |9tLeft), the blue dashed curve shows the modulus of Eq. (IA81 l. If F{(o) is a delta function 5(a) — Olpeak) the spectrum 
obtained by DPT is 



(A9) 


for 0) given by the mesh points {a)jt}. If Olpeak is one of the mesh points, the DPT spectrum has a single peak. However, if Olpeak 
is located in an interval of mesh points, the DPT spectrum has several peaks as shown in Pig. |9jRight). The oscillation of H{co) 
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even gives peaks with negative amplitude. The red rods show how the delta function in the correct spectrum at the position of 
green box appears in finite discrete spectrum in the AC and WK methods. In contrast to the AC method, in the WK method 
the spectrum is squared and the spectrum is given by |F(a))p, hence the spectrum is positive by construction. Although Gibbs 
oscillations are present, they give a width due to the finite time window and do not affect the spectrum much. In Fig.l9landfT0l 
we show comparison of H{o}) and \H{a))\^ respectively. Note that in the WK method, we can obtain the discrete spectrum only 
at the points of even k’s. 




FIG. 9. the spectrum of |/f(<n)| : n = 10000, A = 0.5 : (Left) The true delta-peak is located at k = 0. (Right) The true delta-peak is located at 
k = 0J. 




FIG. 10. the spectrum of ^ \H{a))\'^ : n = 10000, A = 0.5 : (Left)The true delta-peak is located at k = 0. (Right) The true delta-peak is located 
at k = 0.7. 


In order to avoid these apparent negative peaks, a window function, often a Gaussian, is introduced: 



(AlO) 


This treatment replaces by 

a = k = 0,±l,±2,... (All) 

The resulting spectrum |G(®)| is depicted in Fig.fTTI 

The parameter a determines the artificial resolution of the spectrum. Within this resolution, the Gibbs oscillations are smeared 
out, hence we can reproduce the spectrum as discussed in section UTO 
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FIG. 11. The spectrum of |G(a))|: n = 10000, A = 0.5, and a = l. (Left) The true delta-peak is located at L = 0. (Right) The true delta-peak 
is located a.\.k = 0.7 (green line). 

The width of the Gauss window a is known to be taken as 

= £-;► =-21n£, (A12) 

where £ is a number of the order of the smallest number of the computer resolution, say £ = 10^ 


2. Conservation of Intensity 

As mentioned in Sec. |Vl the total amplitude of the spectrum is reproduced correctly by the WK method and the AC method, 
despite the presence of Gibbs oscillations caused by performing DFT. In the subsections that follow, we demonstrate this fact 
analytically for both methods. 


a. Wiener-Khinchin method 

In the case where O) is a continuous variable, the total amplitude of the spectrum is independent of the time T because of the 
following relation; 


1 

T 




(A13) 


In numerical calculations, the Fourier transform is discretized and the Gibbs oscillations may seem to bring some deviation of 
the intensity. But in reality, as shown below, a sum rule such as Eq. ( IA131 l holds even in the discrete case and the correct intensity 
can be calculated without being affected by the Gibbs oscillations, independently of the time T. 

Suppose that there is a peak between two successive (Ok values. Then the DFT and the total amplitude of the spectrum are 
given by 


1 

f 



1 N-\ 

Y" -iCOnkT/N 

k=Q 




2kN'^ 


M 

E 

n=—M 


sin 

sin 



2 


c- ^ ^ c- 27r 

o)„ = — + 5, 0<5<—, 


(A14) 

(A15) 


where 5 is the deviation of the peak’s position from a discretized co, and M should be taken sufficiently large to sum up the 
amplitude around the peak. At the same time we assume that M -^N. Then we can use the following relations: 


sin 


C0„T 

2N 


(OnT 

2 


= sin 


= sin 


/Ttn rSN 




—, ne[-M,M]- 


= (-l)"sin 



(A16) 


Sin 


(A17) 
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Thus the sum is given by 


1 M 
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n=—M 


1 A^-1 

-\(0,ikT/ N ^^ 
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2k 
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V 

2k ^ 

k=0 
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Tin ,T5 

N ^ IN 

1 ■ 2;r3 




1 .2(T8\ ( fTSW 1 

= —^sin —— X TTcosec —— = — 




2 / 


2n' 


(A18) 

(A19) 


where we use the formula of infinite series (n+x)’^ ~ (ncosednx)) ^ ll^ . This result is completely in accord with the 

continuous one Eq. (IA13b . 


b. autocorrelation method 


We can also derive a similar result for the AC method. In the case where o is a continuous variable, the total amplitude of the 
spectrum is independent of the time T because of the following relation: 


\L 


‘At 


dtu = 2n 


/: 


sin{(oT) 


d® = 27t. 


TTCO 


(A20) 


This sum rule holds even in the discrete case as shown below. The basic idea of the derivation is almost the same as in the WK 
method, but it should be noted that d® is half of that in the WK method because the time T is essentially two times larger in this 
case. The discrete version of Eq. ( IA20b is given by 
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E 2Re 
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where 


Using the following relations: 


cos ®„ 
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(OnT 
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and 


the sum is given by 


sm 


(0„T 


= sin{nn + T5) = sin(7r«)cos(7’5) +cos(7r«)sin(7’5) = (—l)”sin(7’5), 


(A21) 

(A22) 

(A23) 

(A24) 

(A25) 

(A26) 
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where we use the formula of infinite series Ln=-oc 


(OnT 
2 

f-D" 

-^ = 2sin(5r) X ncossc(T5) = 2%, 


(A28) 

(A29) 

(A30) 


= TTcosecfTrxl lI^ . 
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Appendix B: Statistical properties of the thermal typical state approach 

The underlying idea of the thermal typical state approach is that it is a good approximation to replace Tr A by (<I>|A|<I>) 
where A = is a D x D Hermitian matrix, the random state |<I>) = ^aW) where the 4,’s are complex-valued Gaussian 
random variables and the set {|a)} can be any complete set of orthonormal states for the Hilbert space of dimension D. The 
demonstration that this approximation is indeed useful requires a proof that by averaging over the we recover the correct 
answer Tr A and that the variance of (<t>|A|<I>) is bounded. 


1. Thermal typical state 

Let us put X = where H ~ and T = which implies X =XK Using Eqs. dlB- •EDl we have 

E[(4>|A|4>)] =2a^TrA = 2c7^Tre^^‘^T, (Bl) 

and 

Var((4>|A|4>)) = 4(j'^Tr 

= 4(7‘*Tr =4-a'^Tr (&-P-^y'^ . (B2) 

In our simulations, we use the state defined by 

\^p) = e-P'^^^\^), (B3) 

to compute estimates of thermal equilibrium expectation values. So we can write With A = e^P"^ and 

B = Y we can derive the upper bound of the variance as following; 

War^ {{^p\Y\^p)) = 16(7^ |Tr {AB)^f=\{{AB)\AB)f 

< {AB,AB){{AB)\ {AB)^) = (Tr {ABYABf = {TrA^B^f , (B4) 

where we used the Schwarz inequality because TrCD = {C,D) defines a scalar product in general. 

As Tr {ABY = Tr > 0, it follows that 

Var ((0/3 \Y\^p)) = 4-a'^Tr (e^P-^Y^ ^ < Tcj'^Tr q-^P-^y^ (B5) 

Tr e“2j3=.^y2 

RSD2((<I>/3|T|<I>/3)) < 2 - (B6) 

(Tr e-P-^T) 


2. Bounds on the variance of the estimate of the partition function 

The inequality Eq. (IB6b is very useful to prove that the statistical error on the estimate of the partition function vanishes 
exponentially with the system size. Specializing to T = 1 and noting that {^p |<I>/ 3 ) > 0 we have 

E[{<Pp\<Pp)]=2a^Tre-P’^ = 2a^Z{l5) = 2ah-P^^P\ (B7) 

Var((0/3|0/3)) =4a'^Tre^2^‘^ = 4a'^Z(2j3) (B8) 

rsd((0/3 10/3)) < Q-Pinm-np))^ (B 9 ) 

where E(j3) denotes the free energy of the system at the inverse temperature j3. As the free energy is an extensive quantity, i.e. it 
is proportional to the number of particles, and E(2j3) — E()3) > 0 for j3 > 0, we have and we obtain 

RSD((0/3 \<pp))< Q-Pi^^^Py^iP)) = (BIO) 

On the other hand, the justification of above estimation may not be obvious for the infinite temperature )3 —> 0. But still the 
convergence of the partition function is guaranteed even in this case because by using 

lime^^^(^) =Trl =D, 

/3-fO 


(Bll) 



21 


we find that 

limRSD{{<Pp\<Pp))<^, (B12) 

in agreement with ll^ . Moreover, by the Schwarz inequality we have in general 

|TrXp = 1(1,< (1,1)(X,X) =DTrX^X, (B13) 


and hence 


RSD{{<^p\<Pp))=RSDmx\<P)) > 
from the definition of RSD given by Eq. (l23t . Thus, it follows that 

limRSD((<I>^|4>^))=^. 


(B14) 


(B15) 


From Eqs. ( IB 101 ) and ( IB15b we can conclude that RSD vanishes exponentially with the system size N. Recall that 

we can choose O’ as we like. From Eq. (IB71 i. it is clear that if we choose o = 1 /V^ we have E[(<I >|3 |<I>| 3 )] = Tr 

that is, the norm of the thermal state is, up to statistical fluctuations which vanish as 1 /D, equal to the partition function. 


3. Approximate estimates 


In this subsection we write Z = e P"^ and to simplify the writing (but without loss of generality), we also choose 0=1/v/2. 
The general idea of the approach is that it suffices to generate one thermal typical state |<I>j3) to find good estimates for (T)eq = 
Tr e-P'^YiTr e- pj^Yhe question is if we can prove that the statistical fluctuations are small in some sense. The problem is the 
following: In the simulation we generate a thermal typical state and compute (fl’/ 3 |T|<t>j 3 )/(<I>p|<I>j 3 ). In the present paper, 
we assume that the denominator and numerator are averaged and then we form the ratio. But if we take the ratio first and then 
we take the average, the estimate is different. In the following we study the variance in this case. Although from Eq. (IB5I) we 
know that variances of these two quantities are bounded from above, we do not yet have a bound on the variance of the ratio of 
them. The purpose of this subsection is to address this point. 

A simple method to estimate averages and variances of the ratio is to make use of the multivariate Taylor expansion for the 
average 


X ^ E[x] Cov[x,y] E[x]Var[y] 
E2[y] E^[y] 


(B16) 


where Cov[x,y] = E[xy] — E[x]E[y] and a similar approximation for the variance 


Var 


LyJ 


Var[x] 2 E[x]Cov[x,y] 


E^[y] 


E3[y] 


E2[x]Var[y] 

E4[y] 


Now we generate a thermal typical state and compute ratios 


(B17) 


R{Y,Z) 


i^pl^p) 


(B18) 


To calculate its average and variance according to Eqs. (IB16b and (IB17b . let us prepare the covariance Cov[('I >|3 |T|'I>| 3 ), (<I >|3 IO 13 )]: 


E[(0^|E|fl>^)(fl>^|fl>^)] = f E[^:^^^;^^]{a\ZY\p){b\Z\q) = TrZ^Y + {TrZY) (Tr Z), 

a,p,b,q=\ 

Cow[{<Pp\Y\<^p),{<^p lO^)] = E[{<Pp\Y\<^p){<^p 10 ^)] - E[{<Pp\Y\<^p)]E[{<Pp = Tr Z^T. 

By using this relation, it follows that 


E 


'(^p\Y\^- 

i^pl^p) _ 


Tr ZT Tr Z^T Tr ZT Var[(<I>j 3 \^p )] 
TrZ (TrZ)^'^TrZ (TrZ)^ 

_ Tr ZT Tr Z^T Tr ZT Tr Z^ 

TrZ (TrZ)^''" TrZ (TrZ)^’ 


(B19) 

(B20) 


(B21) 
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The second and third term in Eq. (IB211 i vanish exponentially with the system size because we see 

TrZ^ 


Tr Z^y Tr ZY Tr Z^ 
& 


(TrZ)^ TrZ (XrZ) 


2 — 


< T 


(TrZ)' 


= T e 


-2j3(f{2j3)-F{i3)) 


where ||y|| = max^, |(v/|y| i//) | is the largest (in absolute value) eigenvalue of T. Hence for large D we have 




TrZy 

TrZ 


— (^)eqj 


as expected. Similarly, for the variance we have 
Var 


'{^plYl^pY 

^ Var[{^p\Y\^p)] ^TrZYTrZ^Y ^ 

^TrZT^ 

^VarKOplOp)] 


(JrZf “ TrZ (TrZ)^ 

[ttZ ) 

(TrZ)^ 


Tr(Zy)2 ^TrZTTrZ^y 


(TrZ)' 

TrZ2 

(Tr Zf 


TrZ (TrZ) 


/TrZy\^ TrZ^ 

V TrZ j (TrZ)^ 


Tr{ZYf ^TrZTTrZ^y 


TrZ2 


TrZ TrZ2 


TrZy 

TrZ 


It is easy to find an upper bound for the terms in the curly brackets. We have 


Tr(Zy)2 

/'TvZyY 

TrZy TrZ^y 

TrZ2 ‘ 

V TrZ ) 

TrZ TrZ2 


<4||Tf, 


and hence we find 


(B22) 


(B23) 


(B24) 


(B25) 


Var 


' {<pp\Y\<P^y 

i^pl^p) 


<4||y||2g-2j3{f(2j3)-F{j3))^ 


showing that the variance vanishes exponentially with the system size N. 


(B26) 
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